################################################
# Analysis for:
#
# Jordan, Soren. 2022. ``Practical Extensions of Regression.'' In 
#  Teaching Undergraduate Political Methodology}, eds. Mitchell Brown, 
#  Shane Nordyke, and Cameron G. Thies. Northampton, MA: Edward Elgar Publishing. 
#  DOI: \verb"10.4337/9781800885479.00028
#
# Required files: None
#
################################################

library(tidyverse)
library(ggpubr)

# Make two x variables, going to arbitrarily rename them M/F and education
set.seed(1)
x1 <- rnorm(100, 5, 2) 		# ``education''
x2 <- rbinom(100, 1, 0.25)	# ``male/female''

y <- x1*rnorm(100, 1, 0.1) + (x2*rnorm(100, 5, 0.25) + rnorm(100, 0, 1)) +
	(x1*x2*rnorm(100, 5, 2) + rnorm(100, 0, 1)) + rnorm(100, 0, 1)

y <- y*1000   # ``salary''

model <- lm(y ~ x1 + x2 + x1:x2)
summary(model)

# Define function to create marginal effects values
meplot <- function(model, var1, var2, ci = .95) {
	alpha <- 1-ci
	z <- qnorm(1-alpha/2)
	beta.hat <- coef(model)
	cov <- vcov(model)
	z0 <- seq(min(model$model[,var2],na.rm=T),max(model$model[,var2],na.rm=T),length.out=1000)
	dy.dx <- beta.hat[var1] + beta.hat[length(beta.hat)]*z0	
	out <- data.frame(x = z0, y = dy.dx)
	out
}

# Store values: this will be top panel
me.data <- meplot(model, var1 = "x2", var2 = "x1")

# Create bottom panel: fitted values/predictions
pred.data <- data.frame(x1 = rep(seq(min(x1), max(x1), length.out = 100), 2),
						x2 = c(rep(0, 100), rep(1, 100)))

# Frame together to make ggplot
gg.pred <- data.frame(x1 = rep(seq(min(x1), max(x1), length.out = 100), 2),
						x2 = c(rep(0, 100), rep(1, 100)),
						y = predict(model, pred.data))

gg.pred$Gender <- gg.pred$x2
gg.pred$Gender[gg.pred$Gender == 1] <- "Male"
gg.pred$Gender[gg.pred$Gender == 0] <- "Female"

pred <- ggplot(data = gg.pred, aes(x = x1, y = y, lty = as.factor(Gender))) +
	geom_line() + theme_bw() +
	labs(x = "Educational Attainment (1-10 Scale)", y = "Predicted Value of Salary", lty = "Gender") + 
	theme(legend.position = c(0.15, 0.75))

me <- ggplot(data = me.data, aes(x = x, y = y)) +
	geom_line() + theme_bw() +
	labs(x = "Educational Attainment (1-10 Scale)", y = "Marginal Effect of Being Male\n(Versus Female) on Salary")

pdf("/.../prme-practical.pdf")
ggarrange(me, pred, nrow = 2)
dev.off()


